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Abstract.  The  one-dimensional  free  energy  model  for  ferroelectric  materials  developed 
by  Smith  et  al.  [28—30]  is  generalized  to  two  dimensions.  The  two-dimensional  free  energy 
potential  proposed  in  this  paper  consists  of  four  energy  wells  that  correspond  to  four 
variants  of  the  material.  The  wells  are  separated  by  four  saddle  points,  representing  the 
barriers  for  90°-switching  processes,  and  a  local  maximum,  across  which  180°-switching 
processes  take  place.  The  free  energy  potential  is  combined  with  evolution  equations  for 
the  variant  fractions  based  on  the  theory  of  thermally  activated  processes.  The  model  is 
compared  to  recent  measurements  on  BaTiC>3  single  crystals  by  Burcsu  et  al.  [8],  and 
predicitions  are  made  concerning  the  response  to  the  application  of  in-plane  multi-axial 
electric  fields  at  various  frequencies  and  loading  directions.  The  kinetics  of  the  90°-  and 
180°-switching  processes  are  discussed  in  detail. 
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1  Introduction 

The  capability  of  piezoelectric  materials  such  as  barium  titanate,  PZT  and  PLZT  ceramics  to  both 
actuate  and  sense  is  due  to  the  non-centrosymmetric  nature  of  the  compounds.  In  actuator  mode,  an 
applied  electric  field  changes  the  ionic  structure  of  the  material,  resulting  in  reversible  or  irreversible 
strains.  The  application  of  a  mechanical  stress  also  alters  the  ionic  configuration  generating  the  voltages 
measured  in  piezoelectric  sensors.  The  two  mechanisms  are  called  the  converse  and  direct  piezoelectric 
effects,  respectively.  The  high  sensitivity  and  repeatibility  of  the  piezoelectric  effects  make  the  mate¬ 
rials  the  current  choice  for  applications  such  as  nanopositioning,  thin-film-based  micro-pumping  and 
sensing.  However,  the  polar  mechanisms  which  provide  piezoelectric  materials  with  the  dual  converse 
and  direct  effects  and  extreme  electromechanical  sensitivity,  also  produce  varying  degrees  of  hysteresis 
and  constitutive  nonlinearities.  Both,  the  hysteresis  and  the  constitutive  nonlinearities  must  be  ac¬ 
commodated  for  high  performance  applications  utilizing  piezoelectric  actuators  and  sensors,  and  for 
this  purpose  an  accurate  and  efficient  model  of  the  material  behavior  has  to  be  developed. 

In  principle,  there  have  been  two  approaches  to  modeling  of  the  nonlinear  behavior  of  piezoelectric 
materials,  one  is  a  macroscopic  phenomenological  approach  and  the  other  one  a  mesoscopic  lattice- 
based  approach.  In  the  former,  a  certain  number  of  macroscopic  internal  variables  are  chosen  to 
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represent  experimental  measurements.  An  early  modeling  attempt  in  this  category  was  made  by  Chen 
[9]  and  Chen  and  Madsen  [10].  Their  work  uses  a  single  scalar  internal  state  variable  that  represents  the 
extent  of  alignment  of  dipoles.  A  detailed  account  of  the  thermodynamic  structure  of  phenomenological 
models  for  ferroelectricity  is  given  by  Bassiouny  et  al.  [4,5]  and  Bassiouny  and  Maugin  [6,7].  Their 
thermodynamic  formulation  uses  the  concept  of  internal  variables  to  capture  a  yield  function  that 
defines  the  domain  of  reversible  behavior  within  the  space  of  applied  fields.  Recently,  Cocks  and 
McMeeking  [13]  have  developed  a  phenomenological  model  that  simulates  the  nonlinear  response  of 
ferroelectrics  to  electrical  and  mechanical  loading.  They  use  the  average  state  of  remnant  strain  and 
remnant  polarization  in  the  polycrystal  as  internal  state  variables.  Kamlah  and  Tsakmakis  [18]  use 
the  remnant  strain  and  remnant  polarization  for  the  constitutive  response  of  ferroelectrics.  However, 
they  decompose  remnant  strain  into  two  parts:  one  is  due  to  remnant  polarization  and  the  other  is  due 
to  ferroelastic  switching.  Their  one-dimensional  model  is  extended  to  allow  multi-axial  loading  [19]. 
Following  the  decomposition  of  Kamlah  and  Tsakmakis  [18],  Kim  and  Kwak  [21]  developed  a  one- 
dimensional  model  for  the  nonlinear  behavior  of  a  piezoelectric  wafer  by  adapting  a  rate-dependent 
viscoplastic  mechanical  model. 

The  second  type  of  approach  is  based  on  a  description  of  the  switching  processes  at  the  crystal  lat¬ 
tice  scale.  In  a  tetragonal  lattice  of  perovskite  crystals,  there  are  six  distinct  types  of  states.  Depending 
on  applied  electric  and  mechanical  loading,  each  of  the  six  states  has  a  different  energy  level,  and  the 
lattices  in  unstable  states  tend  to  change  their  polarization  direction  into  the  polarization  direction  of 
more  stable  states.  By  incorporating  these  features  into  a  switching  model,  Hwang  et  al.  [16]  were  able 
to  reproduce  the  main  effects  found  in  the  macroscopic  response  of  ferroelectric  polycrystals,  namely: 
dielectric  hysteresis,  butterfly  hysteresis  in  strain  versus  electric  field,  and  mechanical  nonlinearity. 
Hwang  et  al.  argued  that  switching  occurs  when  the  work  done  by  local  fields  during  a  given  ferro¬ 
electric  switching  event  exceeds  a  critical  value.  They  used  an  averaging  procedure  to  account  for  the 
variety  of  crystallographic  orientations  and  assumed  uniform  stress  and  electric  fields.  However,  their 
model  showed  a  deficiency  in  a  comparison  with  measurements  of  uniaxial  material  response  [22].  In 
order  to  solve  the  inconsistency,  Hwang  et  al.  [17],  Chen  and  Lynch  [11,12],  Huber  et  al.  [14]  incor¬ 
porated  the  local  interaction  between  adjacent  regions  of  material  with  differing  polarization  states, 
which  was  neglected  in  the  original  work  of  Hwang  et  al.  Kim  and  Jiang  [20]  also  proposed  a  similar 
constitutive  and  finite  element  model  adopting  a  Reuss  type  micromechanics  assumption.  Their  model 
consists  of  Helmholtz  free  energy,  switching  criterion  and  switching  kinetics;  the  rate  of  switching  is 
assumed  to  be  proportional  to  the  thermodynamic  driving  force.  All  of  the  above  models  are  based  on 
rate-independent  or  viscoplasticity-type  switching  evolution  equations.  The  mass  fractions  of  different 
types  of  states  of  the  materials  are  chosen  as  internal  variables,  and  their  evolutions  are  phenomeno¬ 
logically  described  as  a  function  of  the  amount  of  energy  released  during  a  switching  process  or  of  a 
thermodynamic  driving  force.  However,  it  has  not  been  shown  yet  that  these  evolution  equations  are 
capable  of  simulating  the  transient  dynamics  behavior  and  the  closure  of  biased  minor  loops  in  the 
hysteresis  response  of  the  material. 

To  solve  this  problem,  Smith  et  al.  [28—30]  recently  regarded  ferroelectric  switching  as  a  thermally 
activated  process  and  developed  a  one-dimensional  model  based  on  kinetic  switching  equations.  In 
this  approach,  motivated  by  the  shape  memory  alloy  model  developed  by  Muller,  Achenbach,  and 
Seelecke  [1,2,23,27],  a  one-dimensional  Helmholtz  free  energy  potential  consisting  of  two  convex  energy 
wells  and  a  concave  energy  barrier  (spinodal  region)  between  the  wells  is  proposed  as  a  function  of 
polarization.  Depending  on  the  applied  electric  field,  dipoles  jump  from  one  energy  well  to  another 
according  to  a  competition  between  thermal  activation  and  energy  barriers.  Though  their  model  has 
been  successful  in  simulating  the  transient  dynamics  and  minor  loop  behavior  of  the  material,  it 
features  only  two  energy  wells,  and  as  a  result  it  can  not  model  90°-switching  processes. 

This  has  been  improved  by  Sahota  [26] ,  who  extended  the  pure  polarization  model  to  include  a  1-D 
strain  component  and  a  corresponding  90°-variant.  The  approach  allows  to  account  for  the  effect  of 
mechanical  stress  in  one  direction;  however,  Sahota  only  gives  a  qualitative  picture  based  on  purely 
convex  energies  without  spinodal  regimes,  and  no  comparison  with  experimental  data  is  made. 

In  the  present  article,  we  propose  an  extension  of  the  model  by  Smith  et  al.  to  the  two-dimensional 
case.  In  contrast  to  the  full  3-D  case,  this  still  allows  a  graphic  representation  of  the  relevant  energy 
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paraelectric  ferroelectric 


Fig.  1.  Lattice  structure  in  a  paraelectric  and  ferroelectric  PbTiOs  crystal,  90°-  and  180°-switching  processes. 


functions  through  surface  and  contour  plots.  In  order  to  discuss  and  illustrate  the  complexities  arising 
in  multiple  dimensions,  we  therefore  treat  the  2-D  case  in  this  paper  before  addressing  the  development 
of  a  full  3-D  model  at  a  later  stage.  We  also  confine  attention  to  the  case  of  pure  electric  loading  and 
study  the  effects  of  mechanical  stresses  in  a  forthcoming  paper  [3] .  We  discuss  the  additional  complex¬ 
ities  introduced  by  the  multi-dimensional  energy  landscape,  e.g.,  the  determination  of  transformation 
trajectories  due  to  the  field-dependent  movement  of  saddle  points  and  their  impact  on  the  transition 
probabilities  for  the  transformation  from  one  phase  to  the  other.  An  extended  set  of  evolution  equations 
accounting  for  the  increased  number  of  variants  is  presented,  which  is  subsequently  solved  numerically 
for  several  cases  of  electric  loading.  The  model  is  compared  to  polarization  hysteresis  loops  recently 
observed  by  Burcsu  et  al.  [8]  for  BaTiC>3  single  crystals,  and  predictions  are  made  for  nrulti-axial 
electrical  loading  and  various  loading  rates. 


2  Free  Energy  Function 

2. 1  Helmholtz  free  energy  function 

Energy  formulations  for  commonly  employed  ferroelectric  materials  can  be  motivated  by  changes  oc- 
curing  in  the  ionic  structure  during  domain  switchings  or  phase  transitions  in  response  to  applied 
electric  and  stress  fields.  For  illustration,  we  focus  on  PbTiC>3  which  is  isostructural  with  the  mineral 
perovskite  (CaTiCU).  The  crystal  lattice  of  a  perovskite  material,  as  shown  in  Figure  1(a),  exhibits  a 
cubic  configuration  at  temperatures  above  the  Curie  point  T c  (paraelectric  phase)  and  a  tetragonal 
form  below  Tq  (ferroelectric  phase).  Specifically,  a  unit  cell  of  the  material  will  have  a  cubic  structure 
at  temperatures  T>Tc  with  the  Ti+4  ion  located  at  the  center  of  the  lattice,  and  it  will  be  in  one 
of  six  tetragonal  structures  at  temperatures  T<Tc  with  the  Ti+4  ion  biased  along  one  of  the  three 
mutually  orthogonal  crystallographic  directions,  see  for  example  Figure  1(b).  Net  polarization  in  the 
lattice  is  zero  in  the  former  case,  while  it  has  a  finite  value  in  the  latter  case.  The  position  of  Ti+4  is 
not  fixed  but  varies  under  an  electric  field.  Under  an  electric  field  opposite  to  the  polarization  direc¬ 
tion  of  the  lattice,  the  Ti+4  ion  moves  in  the  direction  of  the  applied  field  and,  when  the  electric  field 
exceeds  a  critical  value  (coercive  field),  the  net  polarization  of  the  lattice  is  reversed,  Figure  1(c).  This 
type  of  process  is  called  180°-switching,  while  the  application  of  an  electric  field  perpendicular  to  the 
polarization  direction  of  the  lattice,  leads  to  what  is  called  90°-switching  as  shown  in  Figure  1(d). 

Smith  et  al.  [30]  have  proposed  a  one-dimensional  free  energy  potential  that  models  180°-switching 
processes  with  two  energy  wells  and  one  local  maximum  between  the  two  wells,  as  shown  in  Figure  2.  In 
Figure  2,  the  left  plot  shows  the  one-dimensional  Helmholtz  free  energy  density  ^1(Pi),  and  the  center 
and  right  hand  side  plots  are  the  Gibbs  free  energy  densities  gi(Pi]  Ei)  under  increasing  electric  field. 
The  latter  are  related  to  the  Helmholtz  free  energy  by  the  Legendre  transformation  gi  =  ijj1  —  E\Px. 
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¥i(P1)  =  G1(P1;0) 


g,(p,;e,) 


g,(p,;e,) 


Fig.  2.  One-dimensional  Helmholtz  free  energy  ,01  and  Gibbs  free  energy  g\  for  increasing  field  E\. 


The  application  of  an  electric  field  distorts  the  Gibbs  energy  landscape  and  a  dipole  switch  occurs 
when  the  equilibrium  value  determining  the  Ti+4  position  exceeds  the  unstable  equilibrium  due  to  the 
central  0~2  ion  pairs,  as  shown  in  the  center  and  right  plots  of  Figure  2. 

In  the  present  paper,  the  one-dimensional  free  energy  of  Smith  et  al.  is  generalized  to  a  two- 
dimensional  free  energy  potential  that  consists  of  four  energy  wells,  four  saddle  points  and  one  local 
energy  maximum.  In  the  two-dimensional  energy  model,  the  Helmholtz  free  energy  function  per  unit 
reference  volume  is  given  by 


^(Pi,P2)  =  MPi)  +  m)  +aP?Pi, 

(b(Pi  +  Pit)2  if  Pi  <  -Pi, 

V’i(Pi)  =  {  h  (Pi  -  Pr)(P?/Pi  -  Pr )  if  -Pi  <  Pi  <  Pi, 

[biPi-PR?  if  Pf  <  Pi, 

(b  (P-2  +  Pr)2  if  Pi  <  —Pi, 

MP2)  =  {  h  (Pi  -  Pr)(P2/Pi  -  Pr)  ^  -Pi  <  P2  <  Pi, 

\\r1(P2-PR)2  if  Pi<P2, 


(1) 

(2) 

(3) 


where  ip1  and  are  the  one-dimensional  Helmholtz  free  energy  potentials  in  the  respective  P±-  and 
P2-  directions  that  are  used  by  Smith  et  al.  The  term  Pj  (0  <  P/  <  Pr)  denotes  the  positive  value  of 
polarization,  at  which  the  convexity  of  •01  and  changes  along  the  Pi-  and  P2-  axes,  respectively,  Pr 
stands  for  the  polarization  value  at  the  minima  of  the  four  energy  wells,  and  r]  is  a  material  parameter 
corresponding  to  the  inverse  dielectric  susceptibility  of  the  material.  The  quantity  Pr  is  also  called  the 
remnant  polarization,  and  Pi  is  the  polarization,  at  which  the  material  starts  to  switch  upon  reaching 
the  coercive  electric  field.  The  third  term  in  (1)  is  a  coupling  term  responsible  for  the  four  energy 
wells  to  be  located  along  the  Pi  and  P2  axes.  Surface  and  contour  plots  for  the  Helmholtz  free  energy 
potential  given  by  (1)  are  shown  in  Figures  3(a)  and  3(b),  respectively.  In  Figure  3(b),  one  can  see  four 
energy  wells  denoted  by  Wi,  W2,  W3  and  W4,  each  of  which  being  located  at  (Pr,  0),  (0,  Pr),(-Pr,0), 
and  (0,  —  Pr)  in  the  counterclockwise  direction,  respectively.  In  the  present  paper,  the  notations  Wi, 
W2,  W3  and  W4  are  used  to  denote  the  energy  minimum  states  in  the  corresponding  wells  as  well  as 
the  wells  themselves.  The  four  saddle  points  are  denoted  as  Si,  S2,  S3  and  S4  in  the  figure  and  are 
located  at  (Pc,  Pc),  (—Pc,  Pc),  ( —Pc ,  —Pc),  an(i  (Pc,  ~Pc),  respectively,  where  Pc  is  defined  by 


PC  =  Vfa/2 a)(PR  -  Pi) /Pi.  (4) 

In  contrast  to  the  1-D  case,  there  is  no  longer  a  unique  location  at  which  the  transformation  from 
one  phase  to  another  occurs.  In  fact,  due  to  thermal  fluctuations,  there  is  an  infinite  number  of  paths 
across  the  "ridges"  in  the  energy  landscape  that  the  lattice  elements  can  take  for  the  transformation 
from  one  well  into  another.  A  related  problem  has  also  been  treated  by  Puglisi  and  Truskinovsky 
[24,25]  in  the  context  of  a  system  of  transforming  bi-stable  elements.  They  illustrate  the  various 
transformation  trajectories  of  elements  with  different  barriers  and  construct  Peierls  energy  landscapes 
from  the  projections  of  these  trajectories  across  saddle  points.  To  simplify  the  picture  in  our  case,  we 
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(a)  Surface  plot  (b)  Contour  plot 

Fig.  3.  Two-dimensional  Helmholtz  free  energy  plots. 


will  also  assume  that  the  transformations  always  occur  across  the  saddle  points,  as  these  represent 
the  minimum  energy  barriers.  This  assumption  is  motivated  by  the  observation  that,  even  though 
some  thermal  activation  is  certainly  present  in  the  material,  we  are  dealing  with  a  solid  body,  and  we 
can  assume  it  to  be  relatively  small  compared  to  the  barrier  height  of  the  Gibbs  free  energy  potential. 
Hence,  there  will  be  only  negligible  fluctuation  in  the  transformation  trajectories,  and  the  saddle  points 
represent  the  energy  barriers  for  four  different  types  of  90°-domain  switching  processes.  A  local  energy 
maximum  is  located  at  the  center  of  the  plot,  and  it  is  across  this  maximum  that  the  180°-domain 
switching  processes  occur.  Finally,  it  is  to  be  noted  that  the  two-dimensional  free  energy  given  by 
(l)-(3)  reduces  to  the  one-dimensional  model  of  Smith  et  al.,  only  if  the  electric  field  is  applied  in  the 
±Pi-direction  in  Figure  3  and  no  90°-switching  is  allowed.  In  this  case,  only  180°-domain  switching 
occurs  between  dipoles  in  Wi  and  W3  wells,  and  since  P2  =  0  along  the  Pi  axis,  the  two-dimensional 
free  energy  (1)  reduces  to  ip  (Pi,  P2)  =  ipi(Pi)  +  \t}Pr(Pi  —  Pr )• 


2.2  Gibbs  free  energy  function 

The  last  section  discussed  the  details  of  the  Helmholtz  free  energy  function,  which,  due  to  Ei  =  dip /d Pi, 
is  the  primary  constitutive  quantity  to  be  determined  from  (Ei,  Pi ) -diagrams.  The  free  energy  that 
dictates  the  kinetics  of  the  phase  transformations,  however,  is  the  Gibbs  free  energy,  and  it  is  really 
its  energy  barriers  and  minima  that  are  of  primary  interest. 

In  the  two-dimensional  model,  the  work  of  a  dipole  in  an  electric  field  is  combined  with  the  Helmholtz 
energy  ip  given  by  (1)  to  yield 


g(P1,P2-,E1,E2)  =  iP(PuP2)  -  E^  -  E2P2,  (5) 

where  Ei  and  E2  are  the  components  of  applied  electric  field  in  the  P\  and  P2  directions,  respectively. 
If  an  electric  field  of  magnitude  E  is  applied  at  an  angle  9 e  with  respect  to  the  positive  Pi  axis  in  the 
counterclockwise  direction,  then  Ei  =  E cos Oe  and  E2  =  EsihOe-  The  semicolon  notation  used  in  (5) 
indicates  that  g  is  a  function  of  P,  as  shown  in  Figure  3;  the  electric  field  Et  represents  a  parameter 
shaping  the  energy  landscape  by  changing  the  locations  of  local  minima,  local  maximum  and  saddle 
points. 


3  Basic  Model  Equations 

In  order  to  describe  the  kinetics  of  the  switching  processes,  we  consider  a  representative  volume  element, 
composed  of  mesoscopic  lattice  elements,  each  of  which  see  the  energy  landscape  introduced  above. 


6 


Sang-Joo  Kim  et  al. 


Due  to  thermal  activation,  present  in  every  physical  body  at  non-zero  temperature,  the  polarization 
of  a  lattice  element  fluctuates  about  its  equilibrium  value  in  the  specific  well,  which  characterizes 
its  phase  (±180°,  ±  90°  or  W\  —  W±).  When  an  electric  field  is  applied  in  a  certain  direction,  the 
corresponding  energy  barriers  are  lowered,  and  switching  can  occur  from  one  state  to  another. 

The  average  polarization  of  the  representative  volume  element  in  the  i-direction  is  given  by 

4 

Pi  =  Y,xaP?.  (6) 

a=l 


This  expression  is  the  weighted  sum  of  the  average  polarizations  of  all  lattice  elements  P“  in  the  indi¬ 
vidual  wells,  and  the  weighting  factors  are  the  respective  phase  fractions  xa.  The  average  polarizations 
P"  can  be  computed  from  statistical  arguments,  which  are  based  on  the  probability  to  find  a  particular 
polarization  state  (Pi,P2) 


/j,(P1,P2)  =  C  ex  p 


g(P1,P2]EllE2)  \ 
k bT/Vle  J 


(7) 


This  is  a  typical  Boltzmann  expression  with  a  normalization  constant  C,  ks  and  T  are  Boltzmann's 
constant  and  the  absolute  temperature,  g  is  the  Gibbs  free  energy  density  introduced  above,  and  Vle 
is  the  volume  of  a  lattice  element.  The  average  polarization  in  phase  a  is  then  obtained  from 


pot  _ 


Pi  M  (Pi ,  P2)  dPidP2, 


(8) 


with  the  domain  of  integration  extending  according  to  the  limits  used  in  the  energy  definition  (2)  and 
(3).  In  the  limiting  case  of  vanishing  thermal  activation,  the  probabilities  degenerate  to  Dirac-Delta 
functions,  and  the  average  values  coincide  with  the  polarizations  at  the  minima  of  the  energy  wells. 

The  phase  fractions  xa  in  Eq.  (6)  are  determined  from  a  set  of  evolution  equations,  motivated  by 
the  theory  of  thermally  activated  processes.  The  rate  of  change  of  the  phase  fraction  in  energy  well 
WQ  is  given  by 

4 

xa=  (V/3a  XP  -  Ma/3  x<*)  a  =  1,2, 3, 4,  (9) 

/3=l,/3#a 

which  expresses  the  fact  that  xa  can  change  due  to  a  gain  or  loss  of  lattice  elements  from  or  to 
neighboring  wells.  These  gains  and  losses  are  proportional  to  the  number  of  lattice  elements  in  the 
neighboring  wells,  xp,  and  the  phase  fraction  in  the  well  of  interest,  xa.  The  proportionality  factors 
Hap  are  the  transition  probabilities  for  a  switch  from  W^-well  to  W^-well,  and  they  are  computed 
from  statistical  thermodynamics  as 

Vap  =  ~exP  (fcBr/vLg)  ’  «, /?  =  1,2, 3, 4,  a^/3.  (10) 

They  are  the  product  of  the  probability  of  finding  a  lattice  element  on  top  of  the  barrier  between  Wa- 
and  W^-wells,  and  the  frequency  l/rap  with  which  it  attempts  to  cross  this  barrier.  It  is  well  known 
from  statistical  mechanics  that  these  frequencies  depend  on  the  curvature  of  the  energy  well,  and  if 
this  well  is  not  rotationally  symmetric,  the  frequencies  will  be  anisotropic.  For  the  sake  of  simplicity, 
we  will  not  update  the  frequencies,  when  the  positions  of  the  saddle  points  move,  and  we  will  assume 
t 90  =  r iso  =  r  for  two  different  relaxation  times  T90  and  nso,  representing  different  time  scales  for  the 
two  types  of  switching  processes.  The  term  Agap  :=  g(P°[/3;Ej)  —  g(Pia-,Ej)  is  the  difference  between 
the  energy  on  top  of  the  barrier  at  P“^  and  the  energy  at  the  minimum  of  Wa-well,  P“ .  Note  that 
expression  (10)  is  again  the  low  thermal  activation  limit  of  the  transition  probability,  based  on  the 
same  assumption  that  was  already  introduced  with  the  definition  of  the  saddle  point  as  the  actual 
barrier. 

The  saddle  points,  local  maximum  and  local  minima  are  not  fixed  but  move  in  the  (Pi,  P2)-plane 
depending  on  the  applied  electric  field.  Their  locations  can  be  found  by  solving  the  coupled  nonlinear 
algebraic  equations,  dg/dPi  =  0  and  dg/dP2  =  0,  because  the  gradient  of  the  Gibbs  free  energy  has 
to  be  zero  in  the  P\-  and  P2-directions  at  the  extremum  points.  Usually,  multiple  solutions  exist  in 
the  simultaneous  equations,  and  we  determine  the  type  of  extremum  point  from  the  following  criteria: 
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Fig.  4.  Location  of  extremum  points  of  Gibbs  free  energy  under  increasing  electric  fields  applied  in  the  +P\ -direction. 


If  D(P1,P2)  >  0  and  A(Pi,P2)  >  0 
If  D(P1,P2)  >  0  and  A(P2,  P2)  <  0 
If  D(PuP2)  <  0, 


local  minimum  at  (Pi,P2), 
local  maximum  at  (Pi,P2), 
saddle  point  at  {Pi,P2), 


where  D  and  A  are  discriminant  equations  defined  by 


(11) 


d2G 

d2G 

dP? 

dPxdP2 

d2G 

d2G 

0P1dP2 

OP? 

(12) 


The  evaluation  of  these  equations  shows  that  the  characteristics  of  the  energy  landscape  may  change 
drastically  during  the  course  of  a  loading/unloading  process.  Hence,  we  have  to  carefully  analyze  the 
landscape  in  order  to  define  the  correct  energy  barriers  and  the  corresponding  transition  probabilities. 
The  sequence  displayed  in  Figure  4  illustrates  several  degeneracies  that  may  occur,  and  their  treatment 
is  discussed  in  the  sequel. 
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The  figure  shows  the  movement  of  various  extremal  points  while  the  electric  field  increases  in  the 
+Pi-direction.  In  the  figure,  the  solid  line  represents  dG/dPi  =  0  and  the  dashed  line  dG/dP2  =  0. 
Therefore,  the  points  at  which  the  solid  and  dashed  lines  intersect  are  extremal  points.  Figure  4(a) 
shows  the  various  extremal  points  at  zero  electric  field.  At  E  =  E\,  the  saddle  points  S2  and  S3  move 
vertically  toward  the  Pi-axis  and  the  local  maximum  moves  leftward.  They  meet  at  (—Pc,  0),  and  the 
local  maximum  converts  to  a  saddle  point.  After  that,  the  common  saddle  point  moves  leftward  on 
the  Pi-axis  until  it  meets  the  local  minimum  W3  at  (—P/,0),  as  shown  in  Figure  4(c).  During  this 
period  of  time,  the  energy  barriers  for  the  90°-switching  processes  between  the  W2-  (  or  W4-)  well  and 
the  W3-well  as  well  as  for  all  180° -switching  processes  are  evaluated  at  the  common  saddle  point.  If 
the  electric  field  is  further  increased,  the  local  minimum  W3  and  the  common  saddle  point  merge  at 
the  border  between  the  definition  ranges  of  the  piecewise  energy  function  and  subsequently  vanish.  In 
this  case,  the  values  of  minimum  Gibbs  energy  and  energy  barriers  are  evaluated  at  the  same  location 
(—P/,0).  As  a  result,  the  transition  probabilities  calculated  by  (10)  for  the  switching  processes  from 
the  W3-well  to  the  W/-,  W2-  or  W4-wells  have  maximum  values,  i.e.,  /z31  =  n32  =  /r34  =  1/r.  On  the 
other  hand,  the  saddle  point  Si  moves  up  to  (Pc,  +Pj)  and  merges  with  the  local  minimum  W2  when 
the  field  is  increased,  as  shown  in  Figure  4(d).  After  this,  the  local  minimum  W2  and  the  saddle  point 
Si  do  no  longer  exist,  and  the  transition  probability  for  the  90°-switching  from  the  W2-well  to  the 
Wi-well  is  fi 21  =  1/r,  because  the  values  of  minimum  energy  in  the  W2-well  and  the  energy  barrier 
for  90°-switching  between  W2-well  and  Wi-well  are  evaluated  at  the  same  location  (Pc,+P/).  The 
same  argument  holds  for  switching  between  the  W4-well  and  the  Wi-well.  Comparing  Figures  4(c) 
and  4(d),  we  see  that  the  transition  probabilities  p,32  and  /i34  reach  their  maximum  value  1/r  at  about 
E  =  E2,  but  those  of  /x2i  and  /i41  have  the  same  maximum  value  at  about  E  =  E3.  This  means  that 
90°-switching  in  the  later  period  of  a  process  proceeds  at  a  much  slower  speed  than  the  90°-switching 
in  the  early  period. 

It  is  evident  that  the  determination  of  the  transition  probabilities  in  the  multi-dimensional  case 
requires  careful  tracking  of  the  loading  process.  Once  determined,  the  fia/3  are  then  used  to  assemble 

4 

the  RHS  of  the  ODE  system  (9).  Using  the  phase  fraction  constraint  ]P  xa  =  1,  we  can  finally  reduce 

CK  =  1 

the  four  differential  equations  (9)  to  a  set  of  three  coupled  differential  equations: 

X1  =  — (a*12  +  A* 13  +  A* 14  +  P ll)xl  +  (A*21  —  AUl)2^  +  (A*  31  —  M4l)x3  +  M  41 1 

x2  =  (a1  12  —  ^42)^1  —  (A^l  +  A*  23  +  A*24  +  ^42)^2  +  (1*32  —  A1  42)x3  +  A*42)  (1^) 

x3  =  (Ml3  —  ^43)^1  +  (A* 23  ~  !i4z)x2  ~  (M 31  +  A*32  +  A*34  +  A* 43)^3  +  A*43‘ 

Complemented  by  suitable  initial  conditions,  this  ODE  system  can  be  integrated  numerically  for  a 
prescribed  electric  field  input  to  yield  the  phase  fraction  evolution,  and  together  with  the  algebraic 
equation  (6),  the  time-dependent  polarization  can  be  computed.  Note  that  the  character  of  the  ODEs 

introduces  a  natural  time  scale  through  the  common  relaxation  time  r,  which  allows  to  predict  rate- 

dependent  effects  in  a  natural  way. 


4  Results  and  Discussions 

In  this  section,  the  response  of  the  model  described  in  the  previous  sections  is  calculated  for  several 
different  cases  of  pure  electric  loading.  We  consider  a  single  crystalline  material  with  crystallographic 
direction  coinciding  with  the  global  (Pi,P2)  coordinate  frame.  Therefore,  four  kinds  of  dipoles  exist 
in  the  material  with  their  crystallographic  c  axes  in  the  ±P4-  and  ±P2-axes,  respectively.  First,  we 
simulate  a  quasistatic  polarization  hysteresis  loop  of  a  BaTi03  single  crystal  recently  observed  by 
Burcsu  et  al.  [8],  see  Figure  5.  The  material  parameters  in  (14)  are  chosen  to  fit  the  experimental 
observation. 


r  =  1.3  x  10  2  sec,  Pr  =  0.26  C/m2,  r]  =  7  x  10'  J  ■  m/C2, 
V  =  1.5  x  10^24  m3,  P/  =  0.259  C/m2,  a  =  2  x  108  J  •  to5/C4, 
T  =  298  °K. 


(14) 
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Fig.  5.  Quasistatic  polarization  hysteresis  of  a  BaTiOs  single  crystal,  measurements  by  Burcsu  et  al.  [8],  ( Oe  =  0,  /e  = 
0.05  Hz). 


The  applied  electric  field  is  sinusoidal,  and  its  amplitude  and  frequency  are  Eamp  =  1.1  MV/m  and 
Je  =  0.05  Hz,  respectively.  The  angle  of  the  applied  electric  field  is  denoted  by  6e,  where  0e  is 
measured  from  the  +P\  axis  in  the  counterclockwise  direction.  In  Figure  5,  the  solid  line  represents 
the  result  of  the  simulation  and  the  dashed  line  the  measured  hysteresis  loop.  The  hysteresis  loop 
itself  exhibits  a  very  sharp  transition,  typical  for  a  single  crystal,  and  it  is  rather  thin,  indicative 
of  a  low  coercive  field.  Toward  the  end  of  the  transformation,  the  slope  decreases  smoothly  until 
saturation  is  reached.  We  will  show  this  feature  to  be  a  hint  at  the  potential  micro-mechanism  of  the 
transformation,  which  presents  itself  as  a  two-step  process  with  different  kinetics  due  to  the  presence 
of  90°-  and  180°-switching. 

We  will  first  show  that  the  1-D  two-variant  model  by  Smith  et  al.  is  not  able  to  capture  this 
behavior.  The  prediction  of  the  one-dimensional  model  can  be  obtained  from  the  two-dimensional 
model  by  restraining  all  90°-switching  processes,  that  is,  by  allowing  180°-switching  only.  Identical 
material  parameter  values  are  used  for  the  one-dimensional  calculation,  and  in  Figure  6,  we  compare  the 
response  of  the  two-dimensional  model  (solid  line)  with  that  of  the  one-dimensional  model  (dashed  line). 
While  the  2-D  model  reproduces  the  measured  change  in  slope  toward  the  end  of  the  transformation, 
the  1-D  version  predicts  a  constant  value  for  the  coercive  field.  This  can  be  understood  by  comparing 
the  evolution  of  the  phase  fractions  in  both  cases.  Figure  7(a)  shows  these  evolutions  over  an  electric 
field  range  from  —0.1  MV/m  to  0.7  MV/m  in  the  one-dimensional  case,  and  Figure  7(b)  shows  the 
two-dimensional  case.  Comparing  the  two  figures,  we  see  that  the  spill-out  process  of  lattice  elements, 
initiated  once  the  barrier  shielding  well  W3  is  eliminated,  is  identical  in  both  cases,  and  they  follow 
the  same  curve  in  the  (P,  P)-diagram  initially.  In  the  2-D  case,  however,  these  elements  do  not  all 
jump  into  the  Wi-well  but  into  the  two  lateral  wells  W2  and  W4  as  well.  From  here,  they  transform 
further  onto  Wi,  but  in  a  delayed  manner,  see  again  Figure  7(b).  This  is  due  to  the  circumstance  that 
they  are  trapped  in  wells,  of  which  the  barrier  to  Wi  first  has  to  be  eliminated  by  further  increase 
in  field.  In  fact,  inspection  of  the  phase  fraction  diagram  shows  that  the  combined  percentage  of  W2 
and  W4  fractions  is  as  high  as  70%  at  the  maximum,  and  we  conclude  that  the  more  realistic  2-D 
model  predicts  a  two-step  transformation  through  repeated  90°-switching  processes,  which  is  also  well 
in  agreement  with  the  observation  in  Figure  5. 
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Fig.  6.  Comparison  of  one-  and  two-dimensional  hysteresis  curves  ( 0e  =  0,  / e  =  0.05  Hz). 
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(b)  Two-dimensional  model 


Fig.  7.  Evolution  of  phase  fractions  in  (a)  one-  and  (b)  two-dimensional  model  (Qe  =  0,  f e  =  0.05  Hz) 


Next,  we  discuss  the  response  of  the  model  to  electric  fields  applied  in  three  different  directions 
6e  =  0,  7r/8  and  7r/4.  We  set  the  initial  phase  fractions  to  {x\,X2,xz,x^}\t=o  =  {0,0, 1,0},  and 
therefore,  initially,  all  lattices  in  the  material  are  in  the  W3-well.  The  predicted  polarization  hystereses 
during  the  first  loading  cycle,  with  polarizations  calculated  in  the  direction  of  applied  fields,  are  shown 
in  Figure  8,  and  the  corresponding  plots  of  phase  fraction  evolution  in  Figure  9.  In  Figure  8,  the  solid 
line  corresponds  to  the  electric  field  of  9e  =  0,  the  dashed  one  to  9e  =  7t/8,  and  the  dashed  and  dotted 
one  to  9e  =  tt/4.  The  memory  of  the  initial  state  is  quickly  erased,  and  we  would  observe  steady-state 
hysteresis  loops  after  the  first  cycle  in  all  three  cases.  It  is  shown  in  Figure  8  that  the  value  of  maximum 
polarization  is  the  largest  for  the  electrical  load  at  9e  =  0  and  the  lowest  for  the  electrical  load  at 
9e  =  7t/4.  This  is  due  to  the  circumstance  that  the  maximum  phase  fraction  values  depend  on  the 
angle  of  applied  electric  field  as  shown  in  Figures  9(a)  and  9(c).  At  the  largest  positive  electric  field, 
they  are  approximately  {1,0, 0,0}  and  {0.5, 0.5, 0, 0}  for  9e  =  0  and  7r/4,  respectively.  This  reflects 
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Fig.  8.  Hysteresis  curves  at  three  different  loading  angles  (/e  =  0.05  Hz). 


the  fact  that  the  electric  field,  when  applied  under  increasing  angle  from  0  to  tt/4,  lowers  the  depth 
of  the  W2-well,  as  much  as  that  of  the  Wi-well  at  6e  =  tt/4,  resulting  in  the  equal  probabilities  of 
switching  to  the  Wi-  and  W2-wells.  These  results  are  not  compared  to  measurements  by  Burcsu  et  al., 
who  did  not  study  the  effect  of  different  loading  direction,  and  probably  are  most  closely  related  to 
the  observations  by  Huber  and  Fleck  [15],  who  studied  loading  in  fixed  directions,  but  with  specimens 
cut  under  different  orientations. 

Finally,  in  Figure  10,  we  show  the  effect  of  loading  rate  on  the  polarization  hysteresis.  In  this  case, 
too,  there  are  no  measurements  available  for  comparison  yet,  and  the  results  have  more  of  a  predictive 
character.  We  apply  the  electric  loading  at  three  different  frequencies,  i.e. ,  at  /#  =  0.05  Hz,  0.3  Hz  and 
1  Hz,  and  Figures.  11(a),  (b)  and  (c)  show  the  evolution  of  phase  fractions  at  the  respective  loading 
frequencies  vs.  electric  field  in  the  range  from  E  =  —0.1  MV/m  to  E  —  0.7  MV/m.  In  Figure  10,  it  is 
seen  that  the  hysteresis  loop  gets  wider  with  increasing  loading  frequency.  This  effect  is  small  in  the 
initial  phase  of  the  transformations,  but  it  becomes  very  pronounced  during  the  later  period  of  the 
switching  processes,  when  90°-switchings  occur  from  the  W2-  and  W4-wells  to  the  Wr  or  W3-wells. 
Figure  11  shows  that  the  rate  of  90°-switching  from  W2-  and  W4-wells  to  the  W3-well  gets  slower, 
which  reflects  the  inertia  of  the  transformation.  This  behavior  is  dictated  by  the  ratio  of  the  relaxation 
times  Tigo/90  and  the  rate  of  loading  and  has  to  be  determined  from  experimental  data. 


5  Conclusions 

The  one-dimensional  free  energy  model  for  ferroelectrics  by  Smith  et  al.  is  generalized  to  a  two- 
dimensional  model  that  consists  of  four  energy  wells,  four  saddle  points  and  one  energy  maximum.  In 
the  proposed  model,  polarized  mesoscopic  lattice  elements  jump  from  one  energy  well  to  another  across 
barriers  between  neighboring  wells.  Evolution  equations  for  the  four  relevant  variant  phase  fractions  are 
derived  based  on  ideas  from  the  theory  of  thermally  activated  processes  and  statistical  thermodynamics. 
The  corresponding  transition  probabilities  are  the  product  of  a  Boltzmann  term  giving  the  probability 
to  reach  the  top  of  an  energy  barrier  and  an  attempt  frequency,  which  is  the  inverse  of  a  characteristic 
relaxation  time.  The  effect  of  an  applied  electric  field  on  the  multi-dimensional  energy  landscape  and 
the  barriers  in  particular  is  discussed  in  detail,  along  with  the  algorithm  for  the  determination  of  the 
transition  probabilities. 
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Fig.  9.  Evolution  of  phase  fractions  at  three  different  loading  angles  (/ e  —  0.05  Hz) 


The  prediction  of  the  model  is  compared  with  a  recently  observed  polarization  hysteresis  curve 
for  a  BaTi03  single  crystal.  Using  the  material  parameters  determined  from  the  comparison,  the 
response  of  the  model  to  multiaxial  in-plane  electric  field  loading  at  various  frequencies  is  calculated. 
The  simulations  show  that  the  response  of  the  model  strongly  depends  on  the  direction  of  the  applied 
electric  field.  The  largest  polarization  is  obtained  when  the  electric  field  is  applied  at  9e  =  0,  i.e.,  in 
one  of  the  principal  directions  of  the  free  energy  potential,  and  the  smallest  at  9e  =  7r/4.  The  model 
also  predicts  effects  of  loading  rate,  and  it  is  observed  that  the  hysteresis  loops  become  wider  with 
increasing  loading  frequency,  this  effect  being  particularly  remarkable  during  the  later  stage  of  the 
switching  process.  Inspection  of  the  variant  fraction  evolution  reveals  that  switching  takes  place  as 
a  two-stage  process  with  the  kinetics  being  determined  by  the  energy  landscape  and  the  relaxation 
times. 

In  another  paper  [3] ,  the  implementation  of  mechanical  stresses  and  the  extension  to  polycrystalline 
behavior  is  discussed,  and  a  three-dimensional  version  of  the  model  is  in  development. 
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Fig.  10.  Hysteresis  curves  at  three  different  loading  rates  ( 0E  =  0). 
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Fig.  11.  Evolution  of  phase  fractions  at  three  different  loading  rates  ( 6e  =  0) 
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